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Abstract. A Mathematica based program has been elaborated in order to determine the symmetry 
group of a finite difference equation. 

The package provides functions which enable us to solve the determining equations of the related 
Lie group. 



1 Introduction 

Various works have been carried out on the apphcation of Lie group theory to numer- 
ical analysis. Most of them have been devoted to the building of numerical schemes, 
which preserve the symmetries of the original differential equations. 
Olver [IT] and Kim [10] constructed invariant finite difference equations using the 
concept of the moving frame. Discretization techniques, which preserve some sym- 
metries of the original equations, were studied in [12], [13], [13], [15] and [T6] . 
[12] , [13] , [13 ) [ISj and [IB] used the discrete invariants of the Lie group of the orig- 
inal equation to build invariant finite difference schemes. 

Yanenko [3] and Shokin [2] have provided a Lie group analysis applied to finite dif- 
ference equations, by means of a differential approximation and set down conditions 
under which the differential representation of a finite difference scheme preserves 
the symmetry group of the original differential equation. As the differential approx- 
imation is a differential equation, the Lie group theory can be fully applied to this 
notion. 

The calculation of Lie groups of differential equations with pencil and paper is 
tedious and may induce errors. The size of related equations increases with the 
number of the symmetry variables, and the order of the differential equations. A 
large amount of packages have been created using software programs with symbolic 
manipulations, such as Mathematica, MACSYMA, Maple, REDUCE, AXIOM, MuPAD. 
Schwarz [9] wrote algorithms for REDUCE and AXIOM computer algebra systems, 
Vu and Carminati [6J worked on DESOLVE, a Maple program, Herod [7] and Bau- 
mann [8] developed Mathematica programs. 
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We hereafter describe a new symbolic package, which implements the Lie group 
analysis methods for finite difference equations. The computations are based on 
the theory developed by Yanenko and Shokin. The program has been written for 
Mathematica and provides the symmetry group of a differential representation for a 
given finite difference scheme. The method is based on the Mathematica program of 
Cantwell in pp. 

We presently aim at determining the symmetries lost by the discretization and build- 
ing schemes which preserve those symmetries. 

The method for the investigation of local point transformation groups is set out in 
section 2. The contents of the package is detailed in section 3. Implementation for 
classical numerical schemes is exposed in section 4. 

2 Lie group methods 

2.1 Lie group of differential equations 

Consider a system of /*^-order differential equations: 

= 0, A = l,...,g (2.1) 

Denote by the vector, the components of which are partial derivatives of 

order p, namely, uf^'"''"^ = g^^^^X^^ j = 1, . . . ,n, h, . . . , kp e {1, . . . ,m}. 
Denote by ) the independent variables and u = («!,...,«„) the 

dependent variables. 

The group of local point transformations can be written under the form: 

Gr = {x* = (f)i{x, u, a); u* = ipj{x, u,a), i = 1, ... ,171; j = 1, . . . ,n} (2.2) 

Expand the transformations by means of a Taylor series at the zero value of the 
parameter a^: 



Xi \ Oj(y. 



daa 



+ 0{al), a = l, 

a=0 



Uj = Uj + a„ 



a=0 



+ 0{al), a=l,...,r (2.3) 



The derivatives of (pi and (pj with respect to the parameter are smooth functions, 
called infinitesimals of the group Gf . Denote by 4" and 17" the infinitesimals of Gy. 
In order to find the Lie group transformations of the differential system, it is conve- 
nient to search the infinitesimal operators of Gr'. 

d d 

L„ = (x,ii)— + r7"(x,z/) — , z = l,...,m; j = l,...,n; a = l,...,r (2.4) 

{Lq,, a = 1, . . . , r} represents the set of tangent vectors to the manifold Gr at the 
neutral element a = and is a basis of the Lie-algebra of the infinitesimal operators 

of Gr. 
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The determination of the group transformations is reduced to the determination of 
the infinitesimal functions and r/". 

The knowledge of the enables us to determine the point transformations of the 
group Gr by solving the equations: 

dx* du* 

^=^-(x*,u*), ^ = rjfix*,u*), z=l,...,m; j = a = 1, . . . , r (2.5) 

in conjunction with the initial conditions: 

= x„ u*l^^ = uj (2.6) 



In order to take into account the derivative terms involved in the differential equa- 
tion, the Lie algebra vector field is prolonged: 

cr"'*^'^^^ and (j"'^'^! are given by: 



^ a,{kl...ko-l) Til rso T\t-Ol 

- E , , \ W^' = 2.....' (2.8) 



T^^ko ~t dxidxki . . . dxk,_-, T>Xk„ 



V d duj d 

where: -— = + > 



T>x\^ dxk ^ dxh duj 

Denote by Gr ■* the Lie group of point transformations in the space 
8{x,u,u'^^^\u^^'^'^'^\ . . . ^u^^^^^^'^^'>^ ot thc iudependeut variables, the dependent vari- 
ables and the derivatives of the dependent variables with respect to the independent 
ones. 

Definition 2.1 Consider a subset Vt of the Euclidean space £(a;, w, li(^i)^ ^^C^i-fe)^ ^ ^{ki...ki)^ 

n^{{x, u, u^^'\ u^^'-^'^)-T^{x, u, u^'''\ 1^(^=1 =0, A = 1, . . . , g} (2.9) 

n is an invariant subset of the group Gr^ if all the elements of Gr^ transform any 
point of Q into a point of Q. 



Theorem 2.2 The system of /*'*-order differential equations is invariant under the 
y if: 

L«^\ = 0, a = l,...,r; X = l,...,q (2.10) 



group Gr^ if and only if: 
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2.2 Lie group of differential approximations 

The finite difference scheme, which approximates the differential system fl2.ip . can 
be written as: 

A^{x, u,h,Tu) = 0, X=l,...,q (2.11) 

where h = {hi, h2, ■ ■ ■ , hm) denotes the space step vector, and T = (Ti, T2, . . . , Tm) 
the shift-operator along the axis of the independent variables, defined by: 

(2.12) 

Definition 2.3 The differential equation: 

V\x,u,u'^'''\...,u'^'''-''^'^) = J'^{x,u,u^^^\...,u^^^-^'-^) 

s m 

+J]J](/..)'^7^,^(x,«,«(^^),...,«(^-■'='^.^ 

13=1 i=l 

X= l,...,q; I' = max(x,i)l'\,i (2-13) 

is called the ll^-order differential approximation of the finite difference scheme (12.111) . 
In the specific case s = 1, the above equation is called the first differential approxi- 
mation. 

Denote by G[. a group of transformations in the space S{x, u, h): 

G'j. = {x* = (f}i{x, u, a); u* = ipj{x, u, a); h* = ipi{x,u, h, a), i = 1, . . . ,m; j = 1, . . . , n} 



by Lq,' the basis infinitesimal operator of G'/. 



(2.14) 



d 

Lj = Lc, + dx^u^h) — , a = l,...,r (2.15) 



where 



dan 



, a = l,...,r (2.16) 

a=0 

and by Gr a group of transformation in the space S{x, u, h, u^'^^\ . . . , u^^^'"^^'^). 
The /'^'^-prolongation operator of GJ,, La ^ can be written as: 



£<;■) = w+j:^:,.;*-'.'-^ (2.17) 

i=i p=l ou^ 

Theorem 2.4 The differential approximation (12.131) is invariant under the group 
Gr ^ if and only if 

L(,'')pV(x,m,m('=^),...,u('=^-^'')) =0, a = l,...,r; A = l,...,g (2.18) 



or 



s m 

L«.F^ + Li'')(5^5^(/^,)'''7^,^)]|^^_^ = 0, a = l,...,r; A = l,...,g (2.19) 



P=l i=l 
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Equation f l2.19p leads to a linear overdetermined system of partial differential equa- 
tion, with respect to the infinitesimal functions, called the determining equations of 
the Lie group of the differential approximation (12.131) . Our program determines the 
unknown infinitesimal functions. 



3 Computation Methods 

The program is restricted to partial differential approximations of any order involv- 
ing the unknown scalar function u, the independent variables x and t, the viscosity 
z/ and the step size variables h and r. 



3.1 Calculation of the differential approximation 

First, our program calculates the differential approximation from the knowledge of 
the considered finite difference equation and the approximation error. 
The discrete approximations of the dependent variable, involved in the finite differ- 
ence equation, are expanded at a given order by means of their Taylor series. The 
substitution of these Taylor series expansions into the finite difference scheme pro- 
vides the F-form of the differential approximation, which contains derivatives with 
respect to t, x and mixed derivatives with respect to x and t. The F-form does not 
allow to have stability informations and does not yield the correct order of accuracy 
for all the numerical schemes. That's why we have found essential to determine 
directly the Il-form of the differential approximation, which is obtained by replacing 
the partial derivatives with respect to t and mixed derivatives with respect to x and 
t, involved in the F-form, by partial derivatives with respect to x, using the original 
differential equation. 

The differential approximation is written under the form of an analytic function with 
respect to the independent variables, the dependent variable, the step size variables, 
the viscosity and the partial derivarives of the dependent variable. The dependent 
variable and its partial derivatives are considered as independent variables. 



3.2 Estimation of the determining equations 

Consider the infinitesimal functions , r/", and x". is the infinitesimal related 
to the viscosity: 

X" = T^ ,a = l,...,r (3.1) 
where the viscosity transforms as follows u* = 7r{x, u, u, h, a). 

The remaining infinitesimals of the prolongation operator of the considered symme- 
try group are generated according to the formulae (12.81) . 

The invariance condition of Theorem 12.21 provides a partial differential equation 
involving the unknown infinitesimal functions and products of the partial derivatives 
of the dependent variables. 

Equation (I2.19p is solved as an algebraic equation with respect to the partial deriva- 
tives of the dependent variables, handled as independent variables. Denote by w 
the vector, the components of which are these variables. Since the whole equation 
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holds for all the w components, each coefficient in front of the products of the w 
components has to be zero. This leads to a linear overdetermined system of partial 
differential equation, with respect to the infinitesimal functions, called the determin- 
ing equations of the Lie group of the differential system fl2.13p . The overdetermined 
system is simplified by eliminating the redundancies. This step of the calculation re- 
quires the intervention of the user. The resolution of these equations yields explicitly 
the expression of , r]f, , x", « = 1, • • • , z = 1, . . . , m, j = 1, . . . , n. 

3.3 Solving the determining equations 

The techniques used to solve the determining equations come from [T]. The unknown 
infinitesimal functions are expanded by means of a power series expression with 
respect to the symmetry variables x, t, u, i/, h, r. The polynomial expressions 
are substituted into the determining equations. Solving the determining equations 
amounts to finding the solutions of an algebraic equation. Those latter techniques 
enable us to find in most cases the exact expression of the infinitesimals, when 
the power series is truncated, i.e. when the sought infinitesimals does not contain 
transcendental functions (like exp, cos, sin. In, . . . ). 

3.4 Determination of the symmetry group 

The last part of the program provides the infinitesimal function expression, the Lie 
algebra infinitesimal operators, and the corresponding Lie group transformations. 



4 Examples 

Consider the Burgers equation: 



+ - T^u^x = (4.1^ 

This equation admits the 6-parameter symmetry group: 

* Li = — — , space translation 

ox 

d . 

* L2 = TT, time translation 

ot 

d d <9 , , 

* L3 = X— — h 2t— u—, dilatation 

ox Ot ou 

d d d 

* L4 = xt— — h t^— + i—ut + x) — , projective transformation 

Ox ot ou 

d d 

* L5 = — h — , Galilean transformation 

ox ou 

d d 5 . 

* Lq = —t— + u— — h z/— , dilatation 

Ot ou Of 
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Equation (14. ip can be discretized by means of the finite difference schemes: 
* the FTCS (forward-time and centered-space) scheme: 



T 2h 
* the Lax-Wendroff scheme: 



where: 



r 

+ Y 



2/i 



/i2 



'\2)i+2 2 '^K 2 \2)i 



/l3 



* the Crank-Nicolson scheme: 



r 



I 2 \ 2 ^ \ 2 \ 2 



„2\n+l 



+ 



4h 



2/i2 







The first part of the program provides the differential representation of the schemes: 
* FTCS 

1 



T 







* Lax-Wendroff 

1 



T-2 ]{2 ^1^2 

'^t + 7:(W")^' - ^ ^xx + —g-i + —{U^)xxx - -TTTUxxxx = 

2 D iz iz 



* Crank-Nicolson 



93 . 1 



ut + - i^ti.. + rH- + -{gl + «f72). - -Mxx) + - tt; 



2^ '-^ -^-^ ^^/^ \6 2 12 

where ^1 = - {^) ^ + yuxx, 92 = {-9iu) ^+i^{gi) 93 = {-92U-9j) ^+i^{92) ^ 



-Ur, 
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The next steps of the computation of the symmetry group have been reahzed for 
the above schemes, but only the FTCS scheme is illustrated here. 
The prolonged infinitesimal operator is calculated by means of formula fl2.17p and 
fl2.8p with respect to the infinitesimal functions ^f, 4" and x"- The vector 
the components of which are the derivatives of the dependent variable, treated as 
independent variables, can be written as: 

w = {wt, w2t, w3t, w4:t, wx, wxt, wx2t, wxSt, w2x, w2xt, w2x2t, w3x, w3xt, wAx) 

(4.2) 

The differential representation becomes: 

wt + -u wx — V w2x + — ((7^) + /i^ (-U w3x — —u wAx + -w2x wx) = (4.3) 
2 2 6 12 2 

where g2 = w^ w2x — 2vu w3x + i^^ wAx — Au w2x wx + 2u wx"^. 
The next step provides the determining equations, which are linear partial differen- 
tial equations with respect to the unknown infinitesimal functions. Some of them 
yield informations, which need to be entered by the user: 

er = er(x,t), ^2=^2^, v? = Ux,t)+ug,{x,t), 

C^ = Cnx,t,U,h,T), (^ = C^{x,t,U,h,T), (4.4) 

Then the infinitesimal functions are expanded by means of multivariables polynomial 
expressions: 

6 e-i 

i=0 j=0 



i=Q 



V? = J2Y1 ^'^^'y' + ^J2Y1 dijx'y^, (4.5) 

i=0 j=0 i=0 j=0 

e e~i o-i-j e~i-j-k e-i-j-k-i 



i=0 j=0 k=0 1=0 m=0 

-i—j 6—i—j—k 9—i—j — k~l 



c? = E(E(E( E ( E Wi/'^'ftV")))) 

i=0 j=0 k=0 1=0 m=0 

9—i—j — k 0—i—j—k—l 0—i—j—k—l—m 



E(E(E(E( E ( E 9««.™-V"'AV"v»))))) 

i=0 j=0 fc=0 1=0 m=0 n=0 

The part performing the resolution of the algebraic determining equations gives the 
4-parameter symmetry group represented by: 

= ao + aw X, = bo + (2aio - 5^1)^, < = {gi - aio)u, a < 4 (4.6) 
C" = aioh, = (2aio - 9i)h, = 9ii^ 
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The related 4- dimensional Lie algebra is generated by: 



L' 



L' 



L' 



, (ao = 0, 60 = 1, aio = 0,gi = 0) 



d 

— , (ao = l,bo = 0, aio = 0,gi = 0) 

d_ 

dt 

d d d d d 
X— + 2t— - u— + h— + 2r— , (ao = 0, 60 = 0, Oio = 1, ^1 = 0) 

ox OT Oti oh OT 

d d d d 

+ '^TT + + 2^^^; («o = 0, 60 = 0, aio = 0, 5^1 = 1) 
ox ou oh ov 



(4.7) 



The space and time translations preserve the differential representation of the FTCS 
scheme. 

We recall that (LsjLg) are the vectors of the basis of the Lie algebra Ad of the 
Burgers equation dilatation group. 

{L'^,L'^ represent the vectors of the basis of the Lie algebra A'^ of the differential 

approximation dilatation group. A'^ can be expressed as a direct sum of the Lie 

c) d d d d d 

algebra A^ g spanned by the vectors (La = x— + 2t- — u—,Lr = x— + u— + 2i/— ) 

ox at ou ox ou ou 

and the Lie algebra Ah^^ spanned by the vectors [L'^ = h— + 2r— = ^Jn:^' 

Cf ill CJT CJ ill 

It is natural that if the independent and dependent variables are dilated then the 
step size variables undergo the same transformation. 

La and Lp are elements of the span of the setjLs, Lq}. So A'^ can be represented as 
the direct sum of Ad and Ah^r- The FTCS scheme is invariant under the dilatation 
group, the Lie algebra of which can be written as Ad © ^/i,t- 

The discretization by the FTCS scheme brings about the lost of the Galilean trans- 
formation and the projective transformation. The computation for the other schemes 
shows that all the finite difference schemes admit the same symmetry group. 
The Lie group calculation program has been run under the Mathematica 5.2 version 
for Windows on a PC with a Pentium 4, 2.6 GHz processor. Table [1] shows the in- 



Scheme 


differential approximation order 


Time(s) 


FTCS 


4 


395.2 


Lax-WendrofT 


6 


972.2 


Crank-Nicolson 


6 


906.6 



Table 1: Characteristics 

fluence of the order of the differential approximation on the time of the calculation. 
The increase of differential approximation order yields a sharp rise of the time of 
the calculation. 
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